Inferencia sobre observaciones Normales con JAGS

En este apunte mostramos cómo utilizar JAGS para inferir los valores de posición $\mu$ y de dispersión $\sigma$ de una distribución Normal.

Comenzaremos simulando ciertas observaciones Normales con R a partir de valores conocidos de $\mu$ y $\sigma$. Después utilizaremos esas observaciones en JAGS para evaluar si dicho programa puede recuperar los valores paramétricos que, sabemos, "generaron" las observaciones. Estos dos pasos se conocen como la técnica de "recuperación paramétrica", y se utiliza para evaluar el funcionamiento de algún modelo bajo algún paradigma de inferencia (no necesariamente Bayesiano).

1. Simulación de Observaciones

En la célula siguiente le pedimos a R que genere observaciones provenientes de una distribución Normal con parámetros $\mu=53$ y $\sigma=10$. Noten que en esta simulación el parámetro de dispersión que R espera (en las líneas 4 y 6 de la célula) está en unidades de desviación estándar. Después de simular las observaciones y graficar, incluimos algunos descriptivos de nuevo arreglo z; junto con el histograma, estos primeros pasos son muy importantes para comenzar a estudiar las características de cualquier conjunto de datos

La idea de simular observaciones es que ahora conocemos exactamente qué modelo generó los valores en simul_z, y no sólo eso, sino que sabemos cuáles eran los valores "verdaderos" de $\mu$ y de $\sigma$ "responsables" de esas observaciones.

2. Recuperación de Parámetros

Ahora intentaremos recuperar dichos valores paramétricos programando el mismo modelo (Normal) en JAGS y analizando las distribuciones posteriores correspondientes. En concreto, utilizaremos simul_z como si se trataran de mediciones recolectadas en el mundo natural e intentaremos entenderlas utlizando el modelo:

$$\begin{align} \mu&\sim Uniform(-100,100)\\ \sigma&\sim Uniform(0,60)\\ z_j&\overset{iid}{\sim}Gaussian(\mu,desv. est.=\sigma) \end{align}$$

En otras palabras, comenzaremos suponiendo que el centro de la distribución $\mu$ (también llamado "media") puede estar en cualquier posición entre $-100$ y $100$, y que el parámetro $\sigma$ que representa desviación estándar puede ser cualquier valor entre $0$ y $50$.

Para programar este modelo en JAGS seguiremos la misma estructura que en modelos anteriores, aunque poniendo especial atención a la hora de expresar el parámetro de dispersión en el modelo: a diferencia de R, JAGS no parametriza a la distribución Normal en términos de desviación estándar, sino en términos de precisión $\tau$, en donde la relación entre ambas es:

$$\tau=\frac{1}{\sigma^2},$$

o bien:

$$\begin{align} \sigma^2\tau&=\frac{\sigma^2}{\sigma^2}\\\\ \sigma^2\tau&=1\\\\ \sigma^2&=\frac{1}{\tau}\\\\ \sqrt{\sigma^2}&=\sqrt{\frac{1}{\tau}}\\\\ \sigma&=\frac{\sqrt{1}}{\sqrt{\tau}}\\\\ \sigma&=\frac{1}{\sqrt{\tau}} \end{align}$$

es decir, en un modelo normal lo que se conoce como precisión es igual al recíproco de la varianza, o bien al recíproco de la desviación estándar al cuadrado. Las tres etiquetas, o los tres parámetros ($\tau$, $\sigma$, y $\sigma^2$) sirven para medir la dispersión de la curva, aunque cada uno opera en escalas diferentes. Lo importante es que la relación entre ellos es precisa y podemos programarla para extraer e interpretar el parámetro de dispersión en la escala deseada:

En otras palabras, para escribir el modelo original en JAGS es necesario (y será necesario siempre que trabajemos con modelos Normales) programar esta versión ligeramente modificada:

$$\begin{align} \mu&\sim Uniform(-100,100)\\ \sigma&\sim Uniform(0,60)\\ \tau&\leftarrow\frac{1}{\sigma^2}\\ z_j&\overset{iid}{\sim}Gaussian(\mu,\tau) \end{align}$$

a. Análisis de la Calidad del Resultado de JAGS

Los estadísticos $\hat{R}$ y n.eff, y la apariencia visual de las cadenas, sugieren convergencia confiable. Podemos pasar a examinar los resultados de JAGS para decidir si el algoritmo logró "adivinar" los valores arbitrarios de $\mu$ y $\sigma$ con los que generamos las observaciones $z_j$.

b. Análisis de la Pertinencia del Modelo

(Observaciones vs. Distribución Postdictiva)

El primer paso para examinar el resultado de JAGS siempre consiste en comparar la distribución postdictiva del modelo contra las mediciones de las variables observables. Si estas distribuciones se parecen, podemos continuar interpretando el resto de parámetros, pero si no se parecen y difieren radicalmente no es conveniente continuar en tanto que los parámetros restantes no tendrán sentido respecto de las observaciones porque el modelo que los relaciona es incorrecto.

La distribución postdictiva del modelo se parece a la distribución de las "observaciones" $z_j$. Por lo tanto, podemos continuar examinando el resto de parámetros del modelo.

c. Interpretación de los Parámetros del Modelo

Comenzaremos examinando $\mu$, el parámetro de posición de la curva normal. Presentamos la misma gráfica con dos "zooms" diferentes para resaltar los detalles de las conclusiones posteriores. Aparte, en la gráfica de la derecha hemos incluido el "valor verdadero" de $\mu$ con el que generamos el arreglo $z_j$ para argumentar que Bayes/JAGS recupera razonablemente bien el valor de $\mu$ objetivo.

Continuando con el parámetro de dispersión en escala de desviación estándar ($\sigma$), las gráficas siguientes permiten comparar la incertidumbre prior y posterior sobre este parámetro. Nuevamente, el "valor verdadero" de $\sigma$ se encuentra en una zona de alta densidad posterior en el resultado de Bayes/JAGS.

También exploraremos el parámetro de dispersión en escala de precisión ($\tau$) para enfatizar las diferencias entre esta escala y la escala de desviación estándar. Noten los valores del eje soporte y compárenlos con los valores de $\sigma$. Aparte, es importante resaltar que el suponer incertidumbre a priori uniforme sobre la desviación estándar no implica la misma forma de la incertidumbre sobre $\tau$. En este sentido, es muy importante tener presente sobre cuál de los dos parámetros expresamos nuestra incertidumbre inicial, y recordar que JAGS siempre devolverá resultados en escala de precisión. Como usuarios, es nuestro trabajo traducir dichos resultados a escala de desviación estándar.

Este ejemplo es otra forma de comprobar que Bayes/JAGS rastrea el "valor verdadero" de precisión: noten que dicho valor puede computarse a partir del valor verdadero que elegimos para la desviación estándar, utilizando la transformación en la línea 7 de la célula siguiente:

Ejercicio

La base de observaciones dataset_01.csv contiene registros de una variable.

  1. Lee la base desde R, haz un histograma de las observaciones, y extrae los descriptivos básicos.

  2. Suponiendo que las observaciones provienen de un modelo Normal, infiere los parámetros $\mu$ y $\sigma$ correspondientes.

    a. Ajusta el modelo de JAGS si es necesario, sobre todo respecto de los límites apropiados para los no-observables.

    b. Evalúa el resultado del algoritmo.

    c. Evalúa el desempeño del modelo, comparando la distribución postdictiva contra las observaciones originales.

    d. Interpreta las distribuciones posteriores sobre $\mu$ y $\sigma$.

Solución

  1. Leyendo la base de datos, graficando, y calculando descriptivos:
  1. Implementando el modelo en Normal
$$\begin{align} \mu&\sim Uniform(\mathbf{0},\mathbf{200})\\ \sigma&\sim Uniform(0,60)\\ \tau&\leftarrow\frac{1}{\sigma^2}\\ x_j&\overset{iid}{\sim}Gaussian(\mu,\tau) \end{align}$$

a. Los cambios respecto del modelo anterior se encuentran en las líneas 8, 9, 12, 13, 15, 23, 24, 32, 37, y 38 de la célula siguiente.

b. Evaluando el resultado del algoritmo:

c. Evaluando el desempeño del modelo:

c. Interpretando la incertidumbre sobre las variables no-observables en el modelo:

Apéndice: Notas sobre los histogramas

A continuación presentamos algunas notas sobre la función hist() en R base, con énfasis sobre cómo controlar el ancho de barras y cómo extraer la información del histograma para graficarla en otros formatos.

Para controlar los ejes y las etiquetas recomendamos examinar los argumentos ann=F y axes=F en las gráficas del apunte, así como las funciones axis() y mtext() en dichas células.